EmptyDroplets (FDR <= 0.1) + scDblFindersetwd("/media/jacopo/Elements/re_align/MM/PRJNA732205/SAMN19314095/SRR14629349/")
# Load the libraries (from Sarah script + biomart)
library(tidyverse) # packages for data wrangling, visualization etc
library(Seurat) # scRNA-Seq analysis package
library(clustree) # plot of clustering tree
library(ggsignif) # Enrich your 'ggplots' with group-wise comparisons
library(clusterProfiler) #The package implements methods to analyze and visualize functional profiles of gene and gene clusters.
library(org.Hs.eg.db) # Human annotation package neede for clusterProfiler
library(ggrepel) # extra geoms for ggplo2
library(patchwork) #multiplots
library(reticulate)
Load and do the QC for the cellranger data
#list.files(".")
dat <- Read10X(data.dir ="./out/counts_filtered/")
dat <- CreateSeuratObject(dat) # Create the seurat object from the 10x data
kb.initial <- dat@assays[["RNA"]]@counts@Dim[[2]]
cat("Initial number of cells:", kb.initial,
"\nNumber of genes:", dat@assays[["RNA"]]@counts@Dim[[1]])
## Initial number of cells: 9769
## Number of genes: 36601
Empty cells were already filtered, check for % mt RNA and death markers:
# first calculate the mitochondrial percentage for each cell
dat$percent_mt <- PercentageFeatureSet(dat, pattern="^MT.")
# make violin plots
mt_rna = 20
max_counts = 12000
# Check some feature-feature relationships
# % mt RNA vs n Counts, n Features vs n Counts
# Check some feature-feature relationships
# % mt RNA vs n Counts, n Features vs n Counts
VlnPlot(dat, features = c("nCount_RNA", "nFeature_RNA", "percent_mt")) + geom_hline(yintercept=mt_rna, linetype = "dotted")
plot1 <- FeatureScatter(dat, feature1 = "nCount_RNA", feature2 = "percent_mt")
plot1 <- plot1 + geom_hline(yintercept=mt_rna, linetype = "dotted")
plot2 <- FeatureScatter(dat, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot2 <- plot2 + geom_vline(xintercept = max_counts, linetype = "dotted")
plot1
plot2
## cells retained by mt RNA content ( 20 %): 7308
## percentage of retained cells: 74.81 %
## cells retained by counts ( 12000 ): 7288
## percentage of retained cells: 74.6 %
Check the distribution of the cells with low counts and control death markers:
min_counts = 200
hist(dat@meta.data$nCount_RNA, breaks = 100, xlab = "Counts")
hist(dat@meta.data$nCount_RNA, breaks = 100, xlab = "Counts", xlim = c(0,5000))
hist(dat@meta.data$nCount_RNA, breaks = 1000, xlab = "Counts", xlim = c(0,1000))
abline(v=min_counts, col="red", lty = 3)
The evident peak of cells with < 200 counts could contain dying
cells.
# Subset the dataset to focus only on those cells with low counts
dat.lowcount <- subset(dat, subset = nCount_RNA < min_counts)
# Get the mean of the counts for each gene and sort them decreasing
meanCounts <- rowMeans(GetAssayData(object = dat.lowcount, slot = 'counts'))
meanCounts <- sort(meanCounts, decreasing = T)
# A boxplot can help to observe the distribution of the means
#boxplot(meanCounts)
# Print the most highly expressed genes
head(meanCounts, 30)
## IGKC IGHG1 IGHV4-61 JCHAIN MT-CO2 IGHG3 RPLP1 MALAT1
## 9.4219653 6.8815029 5.2803468 4.8294798 3.4653179 2.5433526 1.8497110 1.7745665
## B2M IGHG4 MT-ND3 MTRNR2L12 RPL34 RPS28 MT-ND4L RPL39
## 1.3786127 1.3699422 1.1213873 1.0809249 0.9884393 0.9161850 0.9161850 0.9132948
## MT-ATP6 RPL32 RPS12 EEF1A1 RPS8 RPL41 MT-CYB SSR4
## 0.8381503 0.7167630 0.7109827 0.6936416 0.6618497 0.6416185 0.6416185 0.6213873
## RPL18A MT-CO3 MT-ND5 RPL10 RPS19 MT-CO1
## 0.6156069 0.5953757 0.5491329 0.5317919 0.5289017 0.5260116
## cells retained by counts ( 200 ): 6942
## percentage of retained cells: 71.06 %
dir.create("result")
saveRDS(dat, file = "./result/SAMN19314095_clean_QC.Rds")
#Normalize
dat <- NormalizeData(dat)
# Find the first 4000 variabe features
dat <- FindVariableFeatures(dat, selection.method = "vst", nfeatures = 4000)
Set mean expression to 0 and variance across 1 to avoid highly expressed genes drive the forwarding analyses. Since negative expression is meaningless, scaled data are useful only for UMAP and clustering
# scale data, the scaled data are saved in:
# dat[["RNA"]]@scale.data
all.genes <- rownames(dat)
dat <- ScaleData(dat, vars.to.regress = c("percent_mt","nCount_RNA"))
dat <- RunPCA(dat, features = VariableFeatures(object = dat), verbose = F, seed.use = 1)
print(dat[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1
## Positive: IGKC, IGHG1, IGHV4-61, JCHAIN, IGHG3
## Negative: IGLL1, STMN1, TYMS, PCLAF, VPREB1
## PC_ 2
## Positive: MS4A1, CD52, HLA-DQB1, HLA-DPA1, HLA-DRB1
## Negative: IGHG1, IGHG3, TYMS, IGKC, IGLL1
## PC_ 3
## Positive: FOS, GADD45B, JUNB, HIST1H2AC, HIST1H2AE
## Negative: B2M, IGHA2, CD27, IGHV3-23, H2AFZ
## PC_ 4
## Positive: CD27, IGHA2, IGHG2, IGHV3-23, IGHA1
## Negative: IGHV4-61, JCHAIN, IGHV4-34, IGKV1-5, PTMA
## PC_ 5
## Positive: AFF3, TCL1A, TOP2A, CFAP73, RUBCNL
## Negative: HCST, S100A4, AIF1, TYROBP, LST1
UMAP is a graph-based method of clustering. The first step in this process is to construct a KNN graph based on the euclidean distance in PCA space:
dat <- FindNeighbors(dat, dims = 1:20)
The graph now can be used as input for the function
runUMAP()
dat <- RunUMAP(dat, dims = 1:20, seed.use = 1)
DimPlot(dat, reduction = 'umap', seed = 1)
## QC metrics
## markers